% This code: discretize the map

function [places_grid,gridmap,edges,places_gecon ] = connected_country_grid( state_borders,places_sedac,...
                                                                                             relief_etopo,places_gecon, Nquantile, group )
    
gridmap = state_borders;    

%% assign altitude and ruggedness to each cell

% X and Y coordinates in topo file
places_X_topo = cell2mat({relief_etopo.centerX}');
places_Y_topo = cell2mat({relief_etopo.centerY}');

% topo variables
av_alt = cell2mat( {relief_etopo.av_alt} );
sd_alt = cell2mat( {relief_etopo.sd_alt} );
rugged = cell2mat( {relief_etopo.rugged} );

for j=1:length( gridmap )

    % roughness
    in_j = inpolygon( places_X_topo, places_Y_topo, gridmap(j).X,gridmap(j).Y  ); % indicate whether each 0.25 degree cell in Etopo belongs in the cell

    gridmap(j).av_alt = mean( av_alt( in_j ) );  
    gridmap(j).sd_alt = mean( sd_alt( in_j ) );  
    gridmap(j).rugged = mean( rugged( in_j ) ); 
    
end


% X and Y coordinates of each city
places_X = cell2mat({places_sedac.X}');
places_Y = cell2mat({places_sedac.Y}');
pop_mat = cell2mat( {places_sedac.population}' ); % vector with city population

% recompute areas and population
for j=1:length( gridmap )      
    
    [   j length( gridmap )       ]

    % Restrict population data to only that country's data 
    % c_pop = arrayfun(@(x) strcmp(cell2mat({x.country}),  icc2name(gridmap(j).Country) ) , places_sedac);

    % total population
    in_j = inpolygon( places_X, places_Y, gridmap(j).X,gridmap(j).Y  ); % indicator that each city is in cell j
    gridmap(j).population = in_j'*pop_mat;  

    % total area
    gridmap(j).area = sum( areaint( gridmap(j).Y,gridmap(j).X ) );

    % if cell remains, compute population centroid
    [ ~,j_largest ] = max( in_j.*pop_mat );  % index of largest city in cell j
    gridmap( j ).coordinates = [ places_X( j_largest ) places_Y( j_largest ) ];
            
end

% Filter out places with 0 population 
gridmap = gridmap(cell2mat({gridmap.population})>0);
             
%% assign income and country to each cell from G-Econ
% assign quantiles
N_cat_income = 20;
quantiles_income = quantile( cell2mat({places_gecon.gdp}),N_cat_income );
for i=1:length(places_gecon)
    [ places_gecon(i).quantiles ] = sum( places_gecon(i).gdp>quantiles_income )+1;
end

% map cells in gridmap to cells in gecon
for i=1:length(gridmap)
    gridmap(i).Xcorner = gridmap(i).X(1);
    gridmap(i).Ycorner = gridmap(i).Y(1);
end

% define a neighbor as a node that is within X miles 
centroidX = zeros(length(gridmap), 1);
centroidY = zeros(length(gridmap), 1);
for j=1:length(gridmap)
    centroidX(j) = gridmap( j ).coordinates(1);
    centroidY(j) = gridmap( j ).coordinates(2);
end

% or neighboring cells (shared border)
border = zeros(length(gridmap));
struct2poly = @(s) polyshape(s.X,s.Y); % define my struct -> polyshape operation
p = arrayfun(struct2poly,gridmap);
pfat = polybuffer(p,1e-4); % was 1e-2 

for k = 1:length(gridmap)
    for j = (k+1):length(gridmap) % Only need j for ks we haven't checked yet
        [ k j ]
        border(j, k) = area(intersect(pfat(j),pfat(k))) > 0;
        border(k, j) = area(intersect(pfat(j),pfat(k))) > 0;
    end
end

for k = 1:length(gridmap)
    gridmap(k).neighbors = [];
    for j = 1:length(gridmap) 
        [ j k ]
        if (border(j,k) == 1 && group == 1) || ...
                (border(j,k) == 1 && group == 2 && ~(j == 83 && k == 20) && ~(j == 20 && k == 83))
                gridmap(k).neighbors = [gridmap(k).neighbors j];
        end
    end
end

%{
for k = 1:length(gridmap)
    gridmap(k).neighbors = [];
    for j = (k+1):length(gridmap) 
        [ j k ]
        if border(j,k) == 1 && (isempty(gridmap(k).neighbors)==1 ...
                || max(ismember(gridmap(k).neighbors, j))==0 )
            if dist_mat(k, j) < 500
                gridmap(k).neighbors = [gridmap(k).neighbors j];
            end
        end
    end
end
%}

% allocate total population to cells in Gecon
for j=1:length(places_gecon)
    [ j ]
    % add up all population in gridmap corresponding to cells match to cell j in Gecon
    places_gecon(j).population_sedac = sum( cell2mat({places_sedac( inpolygon( places_X, places_Y, places_gecon(j).X,places_gecon(j).Y  ) ).population}) );
       
end

% map cells in gridmap to cells in gecon
for i=1:length(gridmap)

    income = 0;
    gridpoly = polyshape(gridmap(i).Y, gridmap(i).X);
   
    % identify the cell in G-econ where the cell i in gridmap belongs to
    for j=1:length(places_gecon)       
        
        geeconpoly = polyshape(places_gecon(j).Y, places_gecon(j).X);
        intersect_region = intersect(geeconpoly, gridpoly);
        
        if intersect_region.NumRegions>0 && places_gecon(j).population_sedac>0
            
            [j places_gecon(j).gdp]
            
            % get population in the intersection area 
            in_pgon = inpolygon( places_X, places_Y, intersect_region.Vertices(:, 2), intersect_region.Vertices(:, 1)  ); % indicator that each city is in cell j
            pop_share = in_pgon'*pop_mat/places_gecon(j).population_sedac;
            
            assert (pop_share<=1.01);
            income = income + places_gecon(j).gdp*pop_share; 
            
        end
        
    end
    
    % income per capita       
    gridmap(i).income = income;
    gridmap(i).y = gridmap(i).income/gridmap(i).population;  %income per capita
    if isnan(gridmap(i).y)
        gridmap(i).y = 0;
        gridmap(i).y = cast(gridmap(i).y,'single');
    end

end

% verify population from G-econ matches with population from SEDAC

% plot( cell2mat( {places_gecon.population} ), cell2mat( {places_gecon.population_sedac} ),'or' )
% corr( cell2mat( {places_gecon.population} )', cell2mat( {places_gecon.population_sedac} )' )

%{
close
mapshow(country_bounds)
mapshow(gridmap(logical(keep)),'FaceColor','none')
%}

% remove new counter as field

%% construct the edges

% vertical-horizontal grid
gridmap_s = [];
gridmap_t = [];
for j=1:length(gridmap)
    gridmap(j).counter = j; % reset the counter bc of the deleted rows
    gridmap_s = [ gridmap_s;gridmap(j).counter*ones( length( gridmap(j).neighbors ),1 ) ];
    gridmap_t = [ gridmap_t;gridmap(j).neighbors' ];
end
temp = [ gridmap_s gridmap_t];

% keep unique horizontal-vertical edges
for i=1:length( temp )    
    temp(i,:) = [ min( temp(i,:) ) max( temp(i,:) ) ];    
end    
[ unique_edges,~,~ ] = unique( temp,'rows');  % we will use unique_edges to build the vertical-horizontal graph

% idenfity cells on boundary
for i=1:length(gridmap)
   if length( gridmap(i).neighbors )<8
       gridmap(i).boundary = 1;
   else
       gridmap(i).boundary = 0;
   end
end

%% assign quantiles

% population
temp =  cell2mat({ gridmap.population })';
quantiles_pop_grid = quantile( temp( temp>0 ), Nquantile  );  % assign quantiles conditional on positive population in the square
for i=1:length(gridmap)
    if temp(i)>0
        [ gridmap(i).quantiles ] = sum( gridmap(i).population>quantiles_pop_grid )+1;
    elseif temp(i)==0
        [ gridmap(i).quantiles ] = 0;
    end
end

% income
temp =  cell2mat({ gridmap.y })';
quantiles_income_grid = quantile( temp( temp>0 ), Nquantile  );  % assign quantiles conditional on positive population in the square
for i=1:length(gridmap)
    if temp(i)>0
        [ gridmap(i).quantiles_inc ] = sum( gridmap(i).y>quantiles_income_grid )+1;
    elseif temp(i)==0
        [ gridmap(i).quantiles_inc ] = 0;
    end
end   

% assign altitude and ruggedness
temp1 =  cell2mat({ gridmap.av_alt })';
temp2 =  cell2mat({ gridmap.rugged })';
for i=1:length(gridmap)
    [ gridmap(i).quantiles_av_alt ] = sum( gridmap(i).av_alt>quantile( temp1, Nquantile  ) )+1;
    [ gridmap(i).quantiles_rugged ] = sum( gridmap(i).rugged>quantile( temp2, Nquantile  ) )+1;
end

%tabulate(cell2mat({gridmap.quantiles})); % check out frequency of population quantiles

%% create a new places file from gridmap
for j=1:length(gridmap)
    
    places_grid(j).Geometry = 'Point';
    places_grid(j).X = gridmap(j).coordinates(1);
    places_grid(j).Y = gridmap(j).coordinates(2);
    places_grid(j).population = gridmap(j).population;
    places_grid(j).income = gridmap(j).income;
    places_grid(j).income = cast(places_grid(j).income,'single');
    places_grid(j).y = gridmap(j).y;
    %places_grid(j).quantiles = gridmap(j).quantiles;  %pop_quantile
    places_grid(j).neighbors = gridmap(j).neighbors;
    places_grid(j).av_alt = gridmap(j).av_alt;
    places_grid(j).sd_alt = gridmap(j).sd_alt;
    places_grid(j).rugged = gridmap(j).rugged;
    places_grid(j).country = gridmap(j).Country;
    
end

edges.unique_edges = unique_edges;
